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We explore the noiseless Burgers dynamics in the inviscid limit, the so-called "adhesion model" 
in cosmology, in a regime where (almost) all the fluid particles are embedded within point-like 
massive halos. Following previous works, we focus our investigations on a "geometrical" model, 
where the matter evolution within the shock manifold is denned from a geometrical construction. 
This hypothesis is at variance with the assumption that the usual continuity equation holds but, in 
the inviscid limit, both models agree in the regular regions. 

Taking advantage of the formulation of the dynamics of this "geometrical model" in terms of 
Legendre transforms and convex hulls, we study the evolution with time of the distribution of 
matter and the associated partitions of the Lagrangian and Eulerian spaces. We describe how the 
halo mass distribution derives from a triangulation in Lagrangian space, while the dual Voronoi- 
like tessellation in Eulerian space gives the boundaries of empty regions with shock nodes at their 
vertices. 

We then emphasize that this dynamics actually leads to halo fragmentations for space dimensions 
greater or equal to 2 (for the inviscid limit studied in this article). This is most easily seen from the 
properties of the Lagrangian-space triangulation and we illustrate this process in the two-dimensional 
(2D) case. In particular, we explain how point-like halos only merge through three-body collisions 
while two-body collisions always give rise to two new massive shock nodes (in 2D). This generalizes to 
higher dimensions and we briefly illustrate the three-dimensional (3D) case. This leads to a specific 
picture for the continuous formation of massive halos through successive halo fragmentations and 
mergings. 

PACS numbers: 47.55.Kf, 98.80.-k, 98.80.Bp, 98.65.-r, 47.27.Gs 



I. INTRODUCTION 

The Burgers equation, introduced at the end of the 
1930s in pQ, was proposed as an archetype system for 
the dynamics of pressureless gases. It indeed shares with 
the Navier-Stokes equations its quadratic nonlinearity, 
and its symmetry properties (invariance under parity, un- 
der space and time translations). It fails however to re- 
produce the chaotic behaviors one encounters in Navier- 
Stokes systems, as it became clear when an explicit fully 
deterministic solution of this system was found by Hopf 
[5] and Cole [3]. Nevertheless, the Burgers dynamics has 
retained much interest for hydrodynamical studies, as a 
benchmark for approximation schemes [4] and as a sim- 
pler example of intermittency phenomena [5HZ]. More- 
over, the Burgers equation has also appeared in many 
physical contexts, such as the propagation of nonlinear 
acoustic waves in non-dispersive media [8] or the forma- 
tion of large-scale structures in cosmology [HI HH] ; see [TT] 
for a recent review. 

In particular, in the cosmological context it provides 
an approximate description of the nonlinear evolution of 
the large-scale structure of the Universe, and it leads to 
complex structures (similar to the network of filaments 
observed between clusters of galaxies in simulations or in 
the sky) in dimension two or larger jTUl fT^HTU] . In this 
cosmological context (and also in many hydrodynamical 
studies) one is interested in the noiseless inviscid limit, 
where the viscosity is sent to zero and when the initial 
velocity field obeys Gaussian statistics. Then, regular 



fluid elements are simply free moving, until singularities 
(shocks) appear. There, the infinitesimal viscosity pre- 
vents matter flows from crossing each other and gives 
rise to shocks and massive local structures of various co- 
dimensions, of infinitesimal width. This effect mimics the 
formation of gravitationally bound objects and allows to 
reproduce the skeleton of the large scale structure built 
by the exact gravitational dynamics with the same initial 
conditions [T31 E] . Then, this model is often referred to 
as the "adhesion model" , following Ref. [9] . 

Many works have been devoted to the early stages of 
the dynamics (see [17] and references in [11]) of such 
models, where only a fraction of the mass has reached 
caustics, and to the universal exponents associated with 
the formation of these singularities. In this paper our 
interest is more particularly focused on a limit where al- 
most all the fluid particles have reached caustics and are 
actually in point-like massive objects. This is the regime 
of interest for cosmological purposes, where the stochas- 
tic initial velocity field is Gaussian with a power-law (or 
slowly running) power spectrum (analogous to fractional 
Brownian motion in ID). Moreover, we are mostly inter- 
ested in the distribution of matter (i.e. the density field 
generated by the Burgers velocity field, starting with a 
uniform initial density), rather than in the properties of 
the velocity field itself. 

Here we must note that the Burgers equation itself, 
and its Hopf-Cole solution, refers to the equation of mo- 
tion for the velocity field. Since we are interested in the 
evolution of the matter distribution, we need to couple 
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this velocity field to a density field and follow their si- 
multaneous evolution. Note that a significant difference 
with the gravitational dynamics is that the evolution of 
the velocity field is independent from that of the density 
field, which in this sense plays a passive role. The "stan- 
dard" procedure would be to use the usual continuity 
equation for the density field, and next take the inviscid 
limit v — > + . This is a natural and well-defined model 
which has been studied for instance in Refs. [TH] HI] ■ A 
drawback of this procedure is that, analytically, it only 
provides information on the instantaneous Eulerian ve- 
locity field inside shocks; to fully obtain the Lagrangian 
map, numerical integration cannot be avoided. One then 
loses the whole interest of having the Hopf-Cole solution 
at our disposal. 

There exists however an alternative procedure for the 
evolution of the matter distribution which is embedded 
in the Hopf-Cole solution for the velocity field. Indeed, as 
noticed in Refs. [Hl[ini[2U], in the inviscid limit the Hopf- 
Cole solution implicitly provides an "inverse Lagrangian 
map" , x i— > q, which for regular points gives the La- 
grangian coordinate q of the particle that is located at 
position x at time t. Once shocks have formed, this map 
is not regular anymore and cannot be inverted without 
ambiguities (except in the one-dimensional case), since 
the "interior" of shock nodes are not reached when the 
whole Eulerian space is spanned. However, using the fact 
that the "inverse Lagrangian map" , x H ¥ q, can be writ- 
ten in terms of a Legendre transform, it is possible to 
define the "direct map" , qn>x, by Legendre conjugacy. 
In fact, this mapping can also be built as the inviscid 
limit of a specific mapping x f-> q, which holds at finite 
v where all points are regular. In this manner, to quote 
Ref. [2U], one obtains an "analytically convenient" model 
for the density field, since the matter distribution can be 
obtained at any time t through Legendre transforms, or 
equivalent geometrical constructions, without the need 
to explicitly solve a differential equation over all previ- 
ous times. Of course, this implies in turn that the den- 
sity field defined by this "geometrical model" does not 
obey the standard continuity equation. More precisely, 
there appears a specific diffusive-type term in the right 
hand side of the continuity equation. Then, in the invis- 
cid limit one recovers the standard behavior outside of 
shocks, but there remains significant differences within 
the shock manifold. This is particularly manifest when 
one tries to follow the history of accretion/merging of 
mass clusters. As already noticed in Ref. [8] (chapter 6, 
see also note [32]), in the geometrical model, and for di- 
mensions two and higher, mass clusters do not necessarily 
merge when they collide, and the collision can give rise 
to several new clusters (with an exchange of matter and 
a possibly different number of outgoing clusters, as we 
shall explain in this paper). This is clearly qualitatively 
at variance with what the standard model gives rise to. 
In the latter case the resulting mass evolution follows in- 
deed that of a genuine "adhesion model" , although some 
unexpected behaviors are still encountered. For instance 



mass clusters can leave shock nodes and travel along the 
shock manifold [18] fT9]. 

The standard model and the geometrical model clearly 
differ and while the first is based on a somewhat natural 
assumption regarding the continuity equation, the sec- 
ond allows easier numerical and analytical insights. In- 
vestigations of the latter model is further justified from 
numerical works that have shown that the latter provides 
a valuable model in the cosmological context, as it builds 
large-scale structure and even shock mass functions that 
are similar to those obtained in the standard gravita- 
tional dynamics [10] . In any case, it provides a rare ex- 
ample of a non-trivial mass transportation, coupled to 
a dynamical velocity field, which can be integrated and 
where significant analytical results can be obtained. 

In this paper we revisit the "geometrical model" , fo- 
cusing on the late time mass distribution after (almost) 
all the mass has reached the shock manifold and gath- 
ered in point-like objects. The Hopf-Cole solution, and 
its geometrical interpretation, provides a way to describe 
the mass distribution in this limit and it corresponds, in 
Eulerian space, to a "Voronoi-like" tessellation (as it was 
already noticed in |14j and further described in more de- 
tails in [THIHO], and in the review paper [3T] ). We further 
explore the consequences of this construction in the con- 
text of cosmological studies by paying special attention 
to the evolution with time of the dual Lagrangian-space 
triangulation and Eulerian-space tessellation. Thus, we 
describe how the halo masses can be obtained, by means 
of Legendre transforms or geometrical constructions, and 
eventually how they are rearranged as time evolves. In 
particular, this allows us to obtain the peculiar collision 
rules that drive the dynamics and to illustrate with fur- 
ther details the exchanges of matter that can take place 
during these events. In this paper we cover those as- 
pects in a rather qualitative way, as the focus will be on 
the description of these mechanisms and their illustration 
with numerical simulations in the two-dimensional case. 
The outline of the paper is the following. The Burgers 
equation and the Hopf-Cole solution are introduced in 
Sect. In] Its dual description in both Eulerian and La- 



grangian spaces is presented in Sect. |III[ together with 
a precise definition of the inviscid limit associated with 
the "geometrical model" and the resulting late time be- 
havior of the displacement fields and potentials. Finally, 
numerical results are presented in Sect. IV and [V] for re- 
spectively the d = 1 and d = 2 cases, as well as a detailed 
description of the fragmentation-merging processes. 



II. BURGERS DYNAMICS 

A. Equations of motion 

We consider the d-dimensional Burgers equation (with 
d>l), 



dtu + (u • V)u = cAu, 



(1) 
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in the inviscid limit, v — > + , for the velocity field u(x, t), 
and the evolution of the density field p(x, t) generated by 
this dynamics, starting from a uniform density p Q at the 
initial time t = 0. As pointed out in the introduction, in 
this article we study a specific "geometrical model" for 
the distribution of matter, where the density field is not 
coupled to the velocity field ([I]) through the standard 
continuity equation but through the modified equation 
(311 below. We shall discuss this "geometrical model" in 
Sect. |III| and explain in details its relationship with the 



Burgers dynamics ([TJ . However, we must first describe in 
this section the initial conditions that we consider in this 
article and the properties of the Burgers velocity field. 

It is a well known result that if the initial velocity is 
potential, u = — V?/>o, it remains so forever [TT], so that 
the velocity field is fully defined by its potential ip(x,t), 
or by its divergence 0(x, t), through 



u = -V^, 6 = -V u = Aip. 



(2) 



Having in mind the use of the Burgers dynamics as a 
model for the formation of the large-scale structure in 
cosmology, we shall consider for numerical implementa- 
tions the case where ipo is a random Gaussian field, sta- 
tistically homogeneous and isotropic, thus entirely char- 
acterized by its power spectrum (see for instance |22| for 
details). The power spectrum of ipo can be equivalently 
defined by that of the velocity divergence, P$ (k), such 
that, 

(0 O ) = 0, (£ (ki)£ (k 2 )) = fo(ki + k a )Pe (fci), (3) 

where we note Sd the Dirac distribution and 9q the 
Fourier transform of the initial divergence, 



o (x)= / dke ik x o (k). 



(4) 



For simplicity we shall also assume that the power spec- 
trum Pg (k) is a power law of index n, 



Pe {k) 



D 



„ri+3— d 



with — 3 < n < 1. (5) 



As a result, for the range —3 < n < 1 the dynamics is 
self-similar: a rescaling of time is statistically equivalent 
to a rescaling of distances, as 



A>0: *->At, x^A 2/(n+3) x, 



(6) 



so that at a given time t there is only one characteristic 
scale, which we normalize as 



L{t) = (2L>i 2 ) 1/( " +3) - 



(7) 



It marks the transition from the large-scale linear regime 
to the small-scale nonlinear regime. In order to express 
the scaling law ^ it is convenient to introduce the di- 
mensionless scaling variables 



Then, equal-time statistical quantities (such as correla- 
tions or probability distributions) written in terms of 
these variables no longer depend on time and the scale 
X = 1 is the characteristic scale of the system. 

We stress however that many properties that we dis- 
cuss in the following also apply to more general initial 
conditions than ([5]), such as those with a scale depen- 
dent power spectrum with a local slope n = d — 3 + 
d In Pg /d In k that could vary with k but remains in the 
range ] — 3, 1[ at very small and very large scales. This 
is typically what is expected in cosmology. 



B. Hopf-Cole solution 

With the Hopf-Cole transformation [2j [3] , tp(x,t) = 
2flnH(x, t), the nonlinear Burgers equation ([!]) trans- 
forms into a linear heat equation for E(x, t), which leads 
to the solution 



^(x,i) = 2z/ln 



dq 



(47Tl4) d / : 



exp 



Then, in the inviscid limit v 
method gives [TTJ |2"3"] 



^o(q) 

2v 

>) 

+ , a steepest-descent 



ipht, t) — max 
q 



Mq) - 



l x ~qr 

2t 



(10) 



If the maximum in (10 1 is reached at a unique point, 



q(x, f), no shock has formed there and the quantity 
q(x, t) is the Lagrangian coordinate of the particle that 
is located at the Eulerian position x at time t [Tl] [53] 
(hereafter, we note by the letter q the Lagrangian coor- 
dinates, i.e. the initial positions at t = of particles, 
and by the letter x the Eulerian coordinates at any time 
t > 0). Moreover, still in the absence of shocks, this 
particle has kept its initial velocity and we have 



u(x,i) = u [q(x,i)] 



q(x,i) 



t 



(11) 



On the other hand, in case there are several degenerate 
solutions to ( 10 ) a shock has formed at position x and 



the velocity is discontinuous (as seen from Eq.( 11 ), as we 
move from one solution q_ to another one q + when we 
go through x from one side of the shock surface to the 
other side) while the density is infinite. 

The solution ( flO| ) has a nice geometrical interpretation 
in terms of paraboloids [TTJ . Thus, let us consider the 
family of upward paraboloids 'Px.cCq) centered at x and 
of height c, with a curvature radius t, 



^x.c(q) 



2t 



(12) 



Then, moving down "P XiC (q) from c = +oo, where the 
paraboloid is everywhere well above the initial potential 
■00 (q) (this is possible for the initial conditions ^ since 
we have |^o(q)| ~ g( 1_,l )/ 2 5 which grows more slowly 
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FIG. 1: The geometrical interpretation of the Hopf-Cole so- 
lution in terms of parabolas for the d = 1 case. 



than q 2 at large distances), until it touches the surface 
defined by ^o(q)) the abscissa q of this first-contact point 
is the Lagrangian coordinate q(x, t). If first-contact oc- 
curs simultaneously at several points there is a shock at 
the Eulerian location x. One can build in this manner 
the inverse Lagrangian map x t— > q(x,t). We show in 
Fig. [TJ an illustration of this geometrical construction in 
the one-dimensional case. 



III. LAGRANGIAN POTENTIAL AND 
MATTER DENSITY FIELD 



A. Finite viscosity v > 

Alternatively, the Burgers dynamics can be entirely 
re-expressed in terms of Lagrangian-Eulerian mappings 
[TU1 ITT] , with the introduction of a time dependent La- 
grangian potential </?(q, t). As we shall discuss below, this 
also allows a precise description of the dynamics from a 
Lagrangian point of view, a precious perspective when 
one is interested in the evolution of the matter distribu- 
tion. 



Note that the inviscid dynamics ([10| is defined as the 
limit of solution (JoJ at v — > + , which is different from 
setting v — in the right hand side of Eq. ([TJ) (that would 
correspond to the so called Zeldovich approximation [24] 
in a cosmological context, which leads to multi-flow re- 
gions if we first go to Lagrangian coordinates). Thus, 
the inviscid limit is singular in this sense and it is use- 
ful to start with v > to remove possible ambiguities 
that would arise by working with Eq.(lO) from the on- 



set. Therefore, following [TO], but for the case of nonzero 
v, let us introduce the function H(x,t), defined as 



H(x,t) 



(13) 



From Eq.Q it also reads as 

if(x,*)=2^1n / e [x-q-|q| 2 /2+^0(q)]/(2,t) 



(14) 



Next, as in [51120], let us define for any function ^4(q) its 



/dqA(q) e [x-q-|q| 2 /2+tV'o(q)]/(2^) 
-<Ay x ^ = J dqe [ x . q -| q |2/2+tV>o(q)]/(2*t) ' ^ 

seen as the average of A over a probability distribution 
for the random variable q, at fixed x and t for a given 
realization of the potential f/'o- This is possible since the 
weight in the average (151 is positive and normalized to 



unity. In particular, we have from Eq.(14| 



dH 
dx 



(16) 



Note that -< q >- x ,t;V>o ^ s n °t the Lagrangian coordinate 
of the particle located at x at time t. Next, the second 
derivatives of H (x) writes as 



d 2 H 



1 



dxidxj 2vt 



2vt 



H3i?j > <<&>"<<&>-] 

<{q— ■«&>•) fe— ■ (17) 



Therefore, in dimension d, the Hessian matrix of H (x) is 
the covariance matrix of the d variables q\, . . . , , which 
implies that it is positive semi-definite [50]. Moreover, 
it has vanishing eigenvalues if and only if the distri- 
bution of the vector q is degenerate (it occurs for in- 
stance when the initial velocity potential is separable, 

"00 (q) = Yn iPo\9i), and V>o ) = for some P air i^o)- 
For generic potentials ipo (l) > such as those obtained from 
the random Gaussian initial conditions pi), the Hessian 
determinant of H(x) is strictly positive, which implies 
that -ff (x) is strictly convex. 

This latter property is crucial. It indeed allows to ex- 
plicitly invert the function q(x, t) =^q>- Xjt; ^, defined in 
Eq.(16). Thus, let us introduce the Legendre transform, 



ip(q, t), of H(x,t) as 

<p(q, t) = £q[H (x, t)] = max [q • x — H(x, t)] 



(18) 



where we used the standard definition of the Legendre- 
Fenchel conjugate /*(s) of a function /(x), 



/* (s) = £ s [/(x)] = sup[s • x - /(x)] 



(19) 



Because H (x) is strictly convex, H (x, t) is also the Legen- 
dre transform of y(q, t) and the correspondence q O x 
associated with the Legendre transforms is one-to-one. 
For a fixed value of q its corresponding value x(q) is 
that which makes q • x — H (x) maximal so that q and x 
obey the relation, 



q(x, t) 



dH 
dx 



(20) 
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As a result q(x) is given by Eq.( 16 1 whereas x(q) is given 
by 



x(q,t) 



dip 
dq- 



(21) 



Finally, note that <p{q) is strictly convex since the Hes- 
sian matrix (d 2 ip/dqidqj) = (dxi/dqj) is the inverse of 
the matrix (dqi/dxj) — (d 2 H/dxidxj), which is positive 
definite, as seen above from Eq.(Jl7fr. 

We are now in position to compare the evolution with 
time of the mapping q x introduced above with actual 
particle trajectories. First, we can note that from the 
definition of the velocity potential |2]) and Eq.(13) the 
Eulerian velocity field also writes as 



uM = -fa = t 



dH 



x- -<q^x,t;y» 

t 



(22) 



where we used Eq.( 16 1. Note the similarity with Eq.( 11 ) 



obtained for the inviscid case at regular points. However, 
in the viscous case u(x, t) is a priori different from the 
initial velocity uo[-<q>~]. It is interesting to note that 
we also have the identity 

x =-<q^ +t -<u (q)>~, (23) 

which arises from a total derivative with respect to q of 



the integrand in Eq.(15), whence 



u(x,t) =-<u (q)^ 



x,t;l/i 



(24) 



However, this does not imply that the Eulerian velocity 
field is set by the initial velocity of Lagrangian particles 
that are located at x at time t. In particular, we have 
~< u (q) uo(^q>-). Nevertheless, we can already see 
from Eqs.(22)-(24| that in the inviscid limit, v — > + , and 
in regular regions, the "mean" (15) becomes dominated 
by a single point q and the variance -< q 2 > c vanishes, in 



agreement with Eq.( 17). We shall then recover the equa- 
tion of motion of free Lagrangian particles. Of course, 
this does not hold within "shocks", where the "mean" 



(15 



is degenerate and receives contributions from sev- 
points q. 

Second, taking the derivative with respect to time of 
the implicit Eq.(16) at fixed -<q>- yields the system of 
equations 

0=^^ + 1^, (25) 



1< i < d : 



dxidxj dt dxidt' 



with summation over repeated indices (here j). Using 
expression (14) to compute the second term and Eq.(17) 
we obtain 







d 2 H 

dxidxj 



dt 



t 



d 3 H 



dxidxjdxj 



(26) 



and multiplying by the inverse of the Hessian matrix, 
(d 2 H/dxidxj)~ 1 , gives 



dxj 
dt 



dxi 
~dt 



= Ui(x, t) — v 



d 2 H 



d a H 



dxi'dxji J ^ dxjdxkdx^ 
(27) 



where we used Eq.(22). In Eq.(27) the first term simply 
changes notation to rewrite dxjdt as a total derivative 
at fixed -< q y, in a form that is more familiar for a 



Lagrangian point of view. Thus, Eq.(27) shows that the 
curves x(-< q y,t), seen as a function of time at fixed 
^.qy, do not follow particle trajectories since their time- 
derivative differs from the Eulerian velocity field u(x, t). 

What is then the resulting density field? As in [THl [2D] , 
we can take advantage of the mapping q o x defined 
by Eq.(16), that is by the Legendre conjugacy (18), to 



construct the distribution of matter at any time t. That 
is, starting with a uniform initial density po a,t t — 0, we 
can define a density field at any time t by 



p(x, t) = po det 



dq 

fa 



which also reads as 



p(x, t) = po det 



2 H 

dxidx. 



Po det 



= po det 



dx 



d 2 V 
dqidq. 



(28) 



(29) 



In Eq.(28) we used the fact that the determinants are 



positive, as shown by expression (29) and the convex- 



ity of i?(x) and y(q), so that we do not need to take 
the absolute value to compute the Jacobians. Of course, 
this model for the density field conserves the total mass. 



Thus, Eq.(29), built from the Legendre conjugacy asso- 
ciated with Eqs.(16)-(21 ), is the definition of the "geo- 



metrical model" that we study in this article. We use the 
label "geometrical" to refer to the fact that it is based on 
Legendre transforms, which have a geometrical interpre- 
tation in terms of convex hulls as we shall discuss below. 



Moreover, it is clear from Eq.(29) that one can build in 



this manner the matter distribution at any time t with- 
out the need to follow the evolution of the density field 
over all previous times, since the convex functions H(x) 
and <p(q) can be computed from the Hopf-Cole solution 
as in Eqs.pTf and pi}. 



Then, from the equation of motion (27), which gives 
the "trajectories" at constant -<q)~, we obtain the evo- 



lution with time of the density field defined by Eqs.(28) 



(29), associated with this mapping x <hM q >-. Thus, 



conservation of mass yields the standard equation 



dp 
dt 



dx 

dt 



0, 



(30) 



which gives with Eq.(27) 
dp 







dt ' V - {pu) = V dx 



p fa- 



de 



/,■/.- 



dxj 

where we introduced the Hessian matrix, O, of H(x) 

d 2 H 



(31) 



(%) 



dxidxj 



(32) 



We can note that in the one-dimensional case, d — 1, this 
is identical to a simple diffusive term Ap, but in higher 
dimensions this is generically different, see the note |43j . 



() 



As proposed in 



the density field denned by the defined as in the viscous case by Eq.(13), we obtain from 



mapping (28) can be interpreted by a stochastic process. Eq.(lO) the expression 



Indeed, with the help of an adequate "mean-field approx- 
imation" one can obtain Eq. ( 28 1 as the density field gen- 



crated by the motion of a passive tracer that moves along 
the Burgers velocity field, to which a Brownian noise of 
amplitude set by v is added (this actually yields a diffu- 
sive term vA.p). If one does not introduce such a "mean- 
field approximation" , an analysis in terms of backward 
stochastic differential equations yields an inverse problem 
with no explicit solution, which only simplifies to ( 28 ) in 
the inviscid limit v — > + |25j . Thus, it might be more 
convenient to simply define the density field through the 
Jacobian (28), or cquivalcntly by the modified diffusion 
equation (31). This provides an explicit interpretation 



in terms of continuum equations of motion as well as 
a simple solution at any finite v in terms of the Leg- 
endre conjugacy (18). As described in the next section, 
this "geometrical model" has a well-defined inviscid limit, 
which allows us to recover the prescription that was used 
in some previous numerical works [lOj . 



B. Inviscid limit v 



We now consider the inviscid limit v — » + . As we 
recalled in Sect. |II B"| the velocity potential is given by the 
maximum ( 10 1, which shows degenerate maxima at shock 



positions. This also determines the Eulerian velocity field 
u(x, t). Next, we define the density field p(x, t) as the 
limit for v — > + of the density field (29) defined in the 
previous section for finite viscosity. 

As already argued in our introduction and in [TBI I20| . 
the main reason supporting the use of the definition ( 29 1, 



that is Eq.(31|, is that it allows an explicit integration 



of the continuity equation as the density field is given 
by Eq.(29) - which remains valid in the inviscid limit as 
seen below - for any time t. In other approaches one 
should a priori numerically integrate the associated con- 
tinuity equation over time. We also stress that the pecu- 
liar form of Eq. ( 31 ) should not be a serious disadvantage, 



since outside of shocks we recover the standard continu- 
ity equation and "within" shocks the Burgers dynamics 
is usually seen as an effective model. 

First, let us define the "linear" Lagrangian potential 
<pL(q,t) by [TO] 



|qr 

2 



#o(q), 



(33) 



so that in the linear regime the Lagrangian map, q n- x, 
associated with particle trajectories, is given by 



Xi(q, t) = = q + tu (q). 



(34) 



Thus we recover the linear displacement field, Xz, — q = 
tfip (q), which is valid before shocks appear, as seen in 
Eq.(|TT|) above. Next, introducing the function -ff(x, t), 



ff(x,t) 



max 
q 



x q 



^|-+#o(q) =£ x [¥>i(q,t)], 

(35) 

where we recognize the Legendre transform of ipi,. As 
recalled below Eq.(10) this gives the inverse Lagrangian 
map x i V q, q(xjT) being the point where the max- 
imum in Eq.( |10p or (35) is reached. Thus, q and x 
are Legendre-conjugate coordinates. From the previ- 
ous section this is also the inviscid limit of the map- 
ping x o-< q y introduced at finite viscosity, and out- 
side of shocks the direct Lagrangian map, q H> x(q, t), 
corresponds to particle trajectories. This mapping is still 
given by Eqs.(16) and (21), which read as 



q(x, t) 



dH 
dx ' 



x(q, t) 



dip 
<9q' 



(36) 



where the function y(q, t) is still defined by the Legendre 
transform (18). From Eq.(|35|) and elementary properties 



of Legendre transforms this implies that ip(q,t) is also 
the convex hull of the linear potential ipL (q) , 



tp(q,t) = Cq[H(x,t)] = conv(ip L ). 



(37) 



Of course, both functions H (x, t) and </j(q, t) remain con- 
vex in the inviscid limit, in agreement with the fact that 
they are still given by Legendre transforms, but they are 
not necessarily strictly convex (i.e. their Hessian deter- 
minant may vanish at some points). By Legendre duality, 
loss of convexity in one of these functions (i.e. hyperpla- 
nar facets) is associated with singular points (with no 
well-defined gradient) in the other one. This analysis 
shows how the prescription (36), that was already used 



in some numerical computations |10j , corresponds to the 
inviscid limit of the density field ( 28 1 and of the mapping 



(16), (21), introduced at finite v in the previous section 



(see the note 44 ) 



Thus, Eqs.rt35fo and (37), which give the convex func- 



tions H(x,t) and ip(q,t), and Eq.(36) which gives the 
Lagrangian to Eulerian space mapping, define the "geo- 
metrical model" that we study in this paper, in the invis- 
cid limit. This provides a simple model for the evolution 
with time of the distribution of matter, which is coupled 
to the Burgers equation for the velocity field. This model 
has already been proposed and studied in previous works 
[H Uni [50], and the main goal of the previous sections 
was to clearly set out its formulation in terms of Leg- 
endre transformations and convex hulls, for both finite 
and infinitesimal viscosity. This also allowed us to derive 



the modified continuity equation (31) to which it corre- 



sponds. In the following we shall investigate some of the 
properties of this "geometrical model" , focusing on the 
regime where all the matter is contained within shocks, 
which holds as soon as t > for the power-law initial 
conditions ([5]). 

The correspondence q f) x is one-to-one as long as 
Vi(q) is smooth and strictly convex (which implies that 
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.ff(x) obeys the same properties). For generic initial con- 
ditions with an ultraviolet cutoff this holds at early times, 
as shown by Eq.([33|. At late times (or also at small 
scales, for the power-law initial conditions that we con- 
sider in this paper), fluctuations of the initial potential 
ipo can make the linear Lagrangian potential tfL non- 
convex. This corresponds to the formation of shocks 



and the maximum (35) is degenerate. In one dimension, 
d = 1, there are generically two degenerate Lagrangian 
points, g_ < q + , which gives rise to a shock at Eule- 
rian position x s . This yields a discontinuity at x s in the 
slope dH/dx(xf) = q± and a discontinuity in the veloc- 
ity u± = (x s — q±)/t. Moreover, since particles cannot 
cross each other (as seen from the convexity of H (x) and 
<p(q), which implies that q(x) and x{q) are monotonically 
increasing) , all the particles initially located in the inter- 
val g+[ map to the shock position x a . 

We can check that this gives for the Lagrangian po- 
tential (p(q), defined by x(q) — dip/dq as in the second 
Eq.([36|, the expression ip(q) = ip^qj) + x s (q— g_) over 
]<?_, q+[. Therefore, tp(q) is indeed the convex envelope ip 
of tp L [10] . Thus, we can see that in the one-dimensional 
case the knowledge of the inverse Lagrangian mapping, 
x n- q, outside of shocks (whence at their boundaries) is 
sufficient to reconstruct the full direct Lagrangian map 
gi->s. This also means that defining the density field as 
the inviscid limit of Eq.( 31 ) (or equivalently in d — 1 with 
a standard diffusive term Ap) , or using the standard con- 
tinuity equation, give the same results. As noticed above, 
this is no longer the case in higher dimensions: the in- 
verse map, x n- q, is no longer sufficient to define the 
direct map, q i— > x, and one must explicitly define the 
model used for the transportation of matter. 



C. "Effective momentum" conservation 

It is interesting to note that the "effective momentum" 
defined by the velocity field u(x, t) is conserved in a global 
sense in the inviscid limit by the "geometrical model" . 
Indeed, let us consider the total "momentum" over an 
Eulerian volume V x , 



dx ;o(x) u(x). 



(38) 



Here we consider a finite viscosity, v > 0, so that the 
velocity and density fields are regular and the integral 
(381 is well defined. Since u is not the actual velocity 
of Lagrangian particles, as explained in sect. |III A[ we 
refer to the quantity defined in Eq.(38) as an "effective 
momentum". Using the expression (28) of the density 
field as the Jacobian of the mapping x f-> q, and the 
result (22), we obtain 



Po 



dq 



*(q) 



t 



dtp 



(39) 



where V q is the Lagrangian space volume which corre- 
sponds to V x through the mapping x O q. Then, con- 



sidering for instance the direction 1 in Cartesian coordi- 
nates, we can integrate over q\ as 



P _ Po 



dq 2 -dq d 



<p(Qi+, ••) ~ ■■) 



(40) 

where qi±(q 2 , ■■,qd) are the boundaries of the volume V q 
at fixed (<?2, •■, <7d) (and we integrate over the projection 
of the volume V q on the (d — 1) remaining directions). 
Here for simplicity we assumed a convex volume (so that 
there are only two boundaries q\±) but the calculation 
straightforwardly extends to the general case. The ex- 
pression (40) holds for any viscosity v. Then, in the 
inviscid limit, v — > + , the Lagrangian potential is given 
by the explicit expression (37) as the convex hull of ip L . 
If shocks are restricted to a finite domain V x hock , which 
corresponds to the finite-mass Lagrangian domain V q hock , 
we can evaluate the total "effective momentum" P over 
a larger volume V x D V. 



shock 



as 



Pi = -Po / dq 2 ..dq d [ip (qi + , ..) - Vo(?i-, ■•)] , ( 41 ) 



where we used Eq.(33), since ip = <Pl a t the surface of V q 



Next, using the definition of the initial velocity potential 
in Eq.([2|, we recognize in Eq.(41) the 1— component of 
the initial momentum Jdqp u o- Therefore, in the invis- 
cid limit the "effective momentum" is conserved in any 
dimension, even after shocks have formed, provided one 
considers the total momentum over a volume such that 
its boundary is regular (i.e. no shock crosses its bound- 
ary). Thus, there is no "dissipation" of the "effective 
momentum" in the inviscid limit. 

In the general case, the "effective momentum" defined 
as the inviscid limit of the expression ( 38 ) may differ from 



the true momentum Jdxpv, where v is the actual ve- 
locity of Lagrangian particles, if shocks have formed. In- 
deed, in regular regions u — > v in the inviscid limit, but 
along shock lines where u(x, t) becomes discontinuous 
these velocities are not identical. (Moreover, v depends 
on whether one uses the "geometrical model" or a "stan- 
dard model" based on the usual continuity equation, as 



discussed in Sect. VD below 



In one dimension, where there is no ambiguity and all 
prescriptions for the density field match in the inviscid 
limit, conservation of momentum in this limit is a well- 
known property |23j . Moreover, in this case the velocity v 
of shock nodes is simply the mean of the left and right ve- 
locities u(x-) and u(x + ), and the "effective momentum" 
( 38 1 is also equal to the actual momentum of Lagrangian 



particles (both are equal to J dqpoUo). 

However, while in one dimension conservation of mo- 
mentum also holds in a local sense, that is, the momen- 
tum of a shock node is equal to the sum of the initial 
momenta of the particles it contains, in higher dimen- 
sions this is no longer the case as shocks can redistribute 
momentum among themselves. We shall come back to 
the conservation of momentum in sect. |VE| below, when 
we discuss the "late-time" regime where all the matter 
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has been redistributed over shock nodes and it is not 
possible to draw volumes with regular boundaries. 

D. Late time structures 

If one lets the system evolve for a long enough time, 
shocks generically form that lead to the formation of low 
dimension structures. Eventually, finite regions V q in 
q-space can map to a single positions x s , leading to a 
point-like massive objects in the Eulerian density field. 
And finally, one expects that, at late enough time, al- 
most all the fluid particles have reached such objects. 
Conversely, it implies that the x-space is partitioned into 
a set of finite domains, V x , each of them originating from 
infinitely small Lagrangian region, q v , centered over a 
discrete number of points. This gives rise to "voids" in 
Eulerian space, as the infinitesimal mass associated with 
Lagrangian position q„ is spread over the finite volume 
V x [EHS] an d this occurs when the initial conditions show 
significant UV power and ipod) has a local maximum at 
q„ with a discontinuous first-derivative. In such regions, 
the function if (x, t) is affine, with a constant gradient 
dH/dx = q„, over the domain V x - 

This is the regime we investigate here. More precisely, 
we assume that </Sz,(q) shows infinitesimally thin down- 
ward peaks over a set of peak locations q^, so that its 
convex hull is supported by a subset of those points. It 
is only in the close vicinity of those points that the map- 
ping is regular, and thus that the local density is there 
independent of the chosen model for the continuity equa- 
tion. However, almost all the mass is then in the shock 
manifold(s) and the question we want to address is how 
the mass is distributed within those manifolds. For the 
ID case the question can be easily answered as shocks 
cover disconnected point-like regions, each of them being 
associated with a single mass halo. For higher dimension 
cases however, the shocks form a single connected mani- 
fold where objects of different dimension and at different 
locations co-exist. It appears that the way the mass is 
distributed within this manifold inherently depends on 
which prescription is adopted for the inviscid limit of the 
continuity equation. 

Let us then explore in more details the case of the 
geometrical model in such a regime. Basic properties of 
its Legendre transform can easily be derived as it then 
reads (for a fixed time t), 

if(x) = max[x-q-<pi(q)] = max[x • q; - ¥>z(qi)]. (42) 

q 1 

It is clear from this form that if (x) is affine over regions, 
corresponding to a given i. Let us label them by the 
corresponding parameter i. Boundaries between two such 
regions, say i and j, correspond to locations where 

x-qi- Vi(qi) = x-qj - <£>i(qj). (43) 

In 2D this obviously corresponds to a straight line (and 
more generally to hyperplanes) . The x-space is then par- 
titioned into polygons, in which if(x) is affine. This 



defines a tessellation, which a priori resembles that of 
Voronoi (we shall see later how the two are linked). 

These regions are simply "voids" where the matter that 
was in the infinitely small q-region % has spread. Note 
that if(x) is continuous, but has discontinuous deriva- 
tives on the boundaries of those regions, which thus cor- 
respond to shocks. Depending on the space dimension 
there are shocks of various dimensions. For instance in 
2D, segments bounding two adjacent domains are shock 
lines. They still gather an infinitesimally small amount of 
matter, which corresponds to that contained in the ridge 
joining two neighbored i regions in q-space (see (TUJ QT] 
for a depiction of those objects). Most of the mass has 
actually reached the point-like shocks (all of the mass at 
any time t > in the regime that we consider in this 
paper). In Lagrangian space they correspond to regions 
where the convex hull y(q) is of a constant slope. 

An important geometrical result is that these regions 
are triangles and that the triangulation is the one asso- 
ciated with the tessellation we have just defined. Indeed 
each point-like shock is associated with each summit x c 
of the tessellation. Those points are at the intersections 
of 3 domains (in 2D) and there are thus three values, 
ii, %2 and 13, such that x c ■ q^ ^hish) i s constant. It 
can easily be checked that no points inside the triangle 
(qii 12, qa) can be part of the convex hull and the whole 
region within the triangle has reached the point x c . In 
terms of mass distribution this means that the mass func- 
tion is given by the area distribution of triangles of that 
triangulation. This property extends to higher dimen- 
sions. In 3D, the masses correspond to the volume of 
tetrahedra [5]. 

This construction shares some similarities with the 
classical Voronoi tessellation/Delaunay triangulation 
dual construction. Let us remind the reader that a 
Voronoi tessellation is built from a set of seeds. It is the 
surface of a volume partition whose domains are such 
that all points of each domain are closest to a given 
seed. The tessellation we have in the present "geometri- 
cal model" is not a true Voronoi tessellation. Indeed, the 
d — 1-hypersurface defining the boundary between two 

domains i and j, , is equivalently defined by 

\<k - x| 2 - 2tVo(q») = |<fe - x| 2 - 2t^ (q J -). (44) 

It corresponds to the hypersurface of the equidistant 
points of q.; and q., only when V'o(qi) = V'o(q.j)- It is how- 
ever possible to embed our d-space into a d+l-space such 
that the extra coordinate q d+1 of q is yj2t{C$ — V'o(q)) 
(where Co is a large enough constant so that all these 
quantities are real, if we consider a finite volume over 
q). Then the d-dimension hypersurface of the equidis- 
tant points between q^ and <\j (in the d + 1 space) is 
defined as 

|q, - x| 2 + (x d+1 - y/2t(C - Mil)))' = 

\qj - x| 2 + (x d+1 - y/2t(C Q -M<H))) ■ (45) 
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Its intersection with the x d+1 = hypersurface is nothing 
but 5y 1 . It straightforwardly follows that the tessel- 
lation we have constructed is nothing but a plane cut 
through a higher dimensional Voronoi tessellation. Note 
that some properties of the Voronoi tessellation have then 
been lost. In particular, seed points may have no domain 
at all associated with (this is when they are not part of 
the convex hull) and the seed point of a given domain 
does not necessarily sit within it. 

It is interesting to note that standard Voronoi tessel- 
lations have also been used in cosmology to study the 
large-scale structures of the Universe. They provide a 
model of these large-scale structures which can reproduce 
some properties of the observed galaxy distribution 
[28] . On the other hand, they can be used as data analysis 
tools, to measure for instance the velocity field statistics 
[59"] . See for instance [301 IS5 f° r reviews. The facts that 
the Burgers dynamics leads to generalized Voronoi cells 
as described above, see also [21], and that this model 
provides a good description of gravitational clustering at 
large scales in cosmology [HI [TO] [TOJ HI] , provide a fur- 
ther motivation for the use of Voronoi tessellations in 
this context. Note that the Burgers dynamics (more pre- 
cisely its formulation as the "geometrical model" stud- 
ied here) also provides the dynamical evolution of these 
tessellations, as a function of the initial conditions. In 
particular, this makes explicit the connection with the 
gravitational dynamics and the background cosmology. 

The picture we have obtained so far is a simple con- 
sequence of the convex hull construction from which the 
geometrical model derives and describes a snapshot of 
the mass distribution. The question we explore hereafter 
is how this structure evolves with time. 



IV. ONE-DIMENSIONAL DYNAMICS 

As time grows one expects the fluid particles to aggre- 
gate in halos of increasing mass. The net result of the 
time evolution should then be that of merging effects. 
For the ID case, this happens simply by the merging of 
adjacent halos. Let us illustrate the mechanism at play 
more explicitly. 



A. Inverse and direct Lagrangian maps 



From the Hopf-Cole solution (10), or the Legendre 
transform ( 35 1 , we obtain the Eulerian velocity field 



u(x,t) and its potential ip(x,t), as well as the inverse 
Lagrangian map, x i— > q. As noticed above, from stan- 
dard properties of Legendre transforms, and as can be 
directly checked on Eq.(35), the map q{x) is monoton- 
ically increasing. This expresses the fact that particles 
cannot cross, because of the infinitesimal viscosity v. We 
plot for illustration in Fig.[2]the results obtained for q(x) 
for one realization of the initial conditions ^ , in the six 
cases n = 0.5, 0, —0.5, —1.5, —2 and —2.5, at some time t. 



O — 




FIG. 2: (Color online) The inverse Lagrangian map, X n> Q, 
for several initial power-spectrum index n. We display the 
results in terms of the dimensionless scaling variables Q and 
X, so that the properties of the curve Q{X) (its increments 
if n < —1) do not depend on time. 



Note that we actually display the dimensionless scaling 
variables X H> Q, defined in Q, so that the time of the 
output is irrelevant since the statistical properties of the 
curve Q(X) (its increments if n < —1) do not depend on 
time. Moreover, we can check that the typical scale (e.g., 
size of vertical jumps or horizontal plateaus) is indeed of 
order unity. 

We can check that Q(X) is monotonically increasing. 
In agreement with previous works [5] [TO] I3"2"H3"5] , we can 
see that for — 1 < n < 1, which we shall label as the "UV 
class" of initial conditions since Pe (k) shows significant 
power at high k, there are large voids (finite intervals over 
X where Q is constant) and a finite number of shocks per 
unit length, where Q(X) is discontinuous and shows pos- 
itive jumps of order unity. For —3 < n < —1, which 
we shall label as the "IR class" of initial conditions since 
Pg (k) shows significant power at low fc, we can distin- 
guish a proliferation of small jumps [6] [TOJ. Indeed, for 
n = —2 it can be shown that shocks are dense in Eule- 
rian space so that there is an infinite number of shocks 
per unit length (the shock mass function diverges at low 
mass) [7, 36 38J. In addition, there are still large shocks, 
of mass of order unity, associated with vertical jumps of 
order unity. 

The initial conditions used to generate Fig. [2] have the 
same Fourier phase for each velocity mode wo(fc) for all n, 
that is, going from the case n\ to n 2 one only multiplies 
each Mo(fc) by fe(™2- n i)/ 2 to satisfy the power spectrum 
([5]). Moreover, the times are chosen so that a given X 
corresponds to the same x (i.e. the figure is obtained from 
different times with n, such that they all have the same 
scale L(t)). As noticed in [TO], this makes the structures 
seen in the various cases very similar, roughly located at 
the same places. This is clearly apparent in Fig. [2] as we 
can see that the various curves Q(X) roughly superpose 
on each other. (For the IR class this is not always the case 
since structures can move over a large distance because of 
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FIG. 3: (Color online) The direct Lagrangian map, Q i— > X, 
for several initial power-spectrum index n. 



the large power at long wavelengths of the velocity field.) 

Next, to derive the direct Lagrangian map, q x, 
we do not need to use the second Legendre transform 
(37) with Eq.(36). Indeed, as noticed in Sect. IIIB in 



one dimension the map q{x) is sufficient to reconstruct 
x(q) since particles do not cross each other, so that both 
functions x(q) and q{x) are monotonically increasing |10j . 
Therefore, spanning q(x) we obtain x(q) (up to the reso- 
lution of our grid). In other words, rotating Fig. [5] by 90 
degrees we can directly read X(Q). We display our re- 
sults in Fig.[3j for the same initial conditions as in Fig. [2] 
We clearly see the horizontal plateaus associated with 
shocks, that form a Devil's staircase for —3 < n < — 1 
where there is a proliferation of small shocks, in agree- 
ment with Refs. [TUIIBI)] . 

For illustration purposes, even though it is not needed 
to derive the Lagrangian map x(q), we display in Fig. [4] 
the linear Lagrangian potential ^>l{q) an d its convex hull 
ip(q), for the cases n = and n — 2, see also [TO]. We use 
the dimensionless scaling variables Q and <E> = ip/L(t) 2 . 
For the white-noise case, n = 0, where tph{q) has no 
finite first derivative, the convex hull only touches <Pl(q) 
at isolated points, which correspond to the boundaries of 
shocks in Lagrangian space (and the position x s of the 
shock in Eulerian space is given by the constant slope 
of (f in-between these boundaries {q^,q + }). Again, we 
can check that typical scales are of order unity. For the 
Brownian case, n — —2, where fL(q) has no finite second 
derivative, we can see some regions of Lagrangian space 
where the convex hull is significantly below ip^ over a 
large interval, associated with a massive shock, but in 
many places it is hard to distinguish both curves. Indeed, 
since shocks arc dense there are finite-size regions which 
contain an infinite number of contact points (note also 
the different scales between the upper and lower panels) . 




Of 




FIG. 4: (Color online) The linear Lagrangian potential <fiL(q) 
and its convex hull p(q), in terms of the dimensionless scaling 
variables Q and $ = <p/L(t) 2 . We show the cases n = 
(upper panel) and n — —2 (lower panel). 



B. Mergings 

As time grows, shocks merge to form increasingly mas- 
sive point-like objects, so that their typical mass scales 
as poL(t). Moreover, once two particles have coalesced 
into a single shock they remain glued together ever af- 
ter. Again, this can be directly seen from the Hopf-Cole 
solution ( fl0| ). Let us have a closer look on Fig. [I] which 
illustrates the parabolic construction (12). On the left 



side we show the case of a regular point x' , associated 
with a unique Lagrangian coordinate q' , whereas on the 
right side the first-contact parabola V x . c simultaneously 
touches the potential ipo (q) at the two first-contact points 
q_ and g+. This implies the points within ]q^,q + [ can- 
not come into contact with a parabola. There is actu- 
ally a shock at the Eulerian position x, which gathers all 
the mass associated with the Lagrangian interval ]q_ , q+ [. 
Then, as the parabola curvature decreases with time as 
1 jt the points within the interval ]q- , q+ [ will remain un- 
reachable. This implies that this interval cannot be bro- 
ken into separate pieces and the whole interval ]q-,q + [ 
can only belong to a single shock, which can eventually 
grow larger. 

We show in Fig. [5] the particle trajectories obtained 
for one realization in the case n = 0, using the physical 
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FIG. 5: (Color online) The particle trajectories obtained for 
one realization in the white-noise case n = 0. For any time, 
each line correspond to a collection of particles that have 
reached the same location. 



coordinates x and t, see also |3S]- We can see how parti- 
cles are gathered into isolated shocks that merge as time 
grows to build increasingly massive and rare point-like 
masses. This merging process happens when two lines 
meet. At such points two halos merge into a single halo. 
The mass gathered in a single halo then corresponds to 
the volume span by all its progenitors. The average mass 
obviously grows with time (according to the scaling law 
^ for the initial conditions Note that in between 

mergings halos move along straight lines, in agreement 
with Eq.Q. At merging point their velocity changes. 
Its new value corresponds to the average of the mass 
weighted velocities of the halos that merge together [33] . 
Note finally that in between lines there is no matter left. 
These are voids. 

Other spectrum indices in the range — 1 < n < 1 yield 
similar figures (with shocks getting more numerous for 
lower n), while for —3 < n < —1 the particle trajectories 
fill the whole (x, t) plane since shocks are dense in the 
Eulerian space. Note however that their boundary lines 
are not dense in Lagrangian space. This is simply due 
to the fact that although their number density is infinite, 
halos have finite masses. 



V. TWO-DIMENSIONAL DYNAMICS 

We now consider the two-dimensional case, d = 2. The 
Hopf-Cole solution and its associated parabola construc- 
tion can be transposed from the ID case. The merging 
properties of the halos are however much more intricate. 



A. Lagrangian and Eulerian tessellations 

In high dimensions the inverse Lagrangian map, x i— > 
q, is still given by the Hopf-Cole solution ( fl0| ), that 
is, by the Legendre transform ( 35 1-( 36 ) . In the one- 



dimensional case, when shocks have formed, that is when 
x maps to the degenerate pair monotonic- 
ity arguments ensure that conversely the whole interval 
q + [ maps to position x. For higher dimensional cases 
this is no more the case : the image of the function q(x) 
does not span the whole q-space so that the map x M> q 
does not fully determine the inverse function x(q) by it- 
self (and a further prescription must be added). That 
would be the case if the first-contact paraboloid ? x . c , 
introduced in Eq.(fl2j), associated with a massive shock 
node at x, would make contact with the potential ipo(<l) 
over a closed curve C x . But this is generically not so: V x , c 
only touches the potential ?/>o(q) over a few points - from 
one for a regular point to d + 1 for a shock node. As de- 
scribed in Sect. |III[ this missing information is provided 
by the Lagrangian potential (p(q), through Eqs.(36) and 
(37), which take into account that the maps q(x) and 
x(q) arise from the Burgers dynamics, to which we have 
coupled the evolution of the matter distribution defined 
by the "geometrical model" . 

Since we focus in this section on the two-dimensional 
case, for generic initial velocity potential ipo the convex 
hull tp is made of regular parts, ruled surfaces and planar- 
triangles, see Fig. 5 of [TUj. However, for the power-law 
initial conditions (5|, where the initial velocity potential 
obeys the scaling law 



A > : V^o(Aq) = A( 1 - n '/ 2 ^ (q), 



(46) 



so that v?i(q) shows no finite first derivative if — 1 < n < 
1, and no finite second derivative if —3 < n < — 1, there 
are no regular parts and the convex hull is entirely made 
of planar triangles. This corresponds to the late time 
limit described in Sect. 



HID 



a regime in which all the 
mass is located within shock nodes, see also [8]. Note that 
this holds for any dimension, for the initial conditions ([5| , 
the convex hull <^(q) being made of pieces of hyperplancs 
of dimension d, defined by d+1 points. This is illustrated 
on Fig. [4] for the one-dimensional case, on Fig. [6] for the 
two-dimensional cases. 

We show in Fig. [6] the triangulations we obtain in the 
q-space at a given time and for one realization of the 
initial potential tp , for the two cases n = (upper panel) 
and n = — 2 (lower panel). We can see that we recover 
the qualitative properties found in the one-dimensional 
case. In the case n = contact points of the convex 
hull (which are the nodes of the triangulations) are well 
separated and most triangles have an area of order unity, 
which implies that in Eulerian space shock nodes are in 
finite number per unit surface with masses of order unity 
(in the dimensionless scaling units ([8])). In the case n = 
—2 there are regions of the Lagrangian space with an 
infinite number of nodes, with numerous triangles of very 
small area (up to the numerical resolution) , which implies 
that in Eulerian space there is also an infinite number of 
shocks per unit surface (in the mean), that is, the mass 
function of shocks diverges at low masses. We can also 
distinguish some "white" areas in the figure, associated 
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n=0 



n=0 




FIG. 6: (Color online) The partition of the Lagrangian space 
associated with triangular facets of the convex hull y>(q) of 
the linear Lagrangian potential. We plot the triangulation 
in terms of the dimensionless scaling coordinates (Qi,Q2)- 
The Lagrangian potential realization corresponds to the case 
n = (upper panel) and n = —2 (lower panel). 



with triangles of area of order unity, that is, shock nodes 
with a mass of order unity. 

The Lagrangian space triangulation shown in Fig. [6j 
which arises from the triangular facets of the convex hull 
</?(q), is associated with a dual partition of the Eulerian 
space, which arises from the planar facets of the Legen- 
dre transform H(x). Indeed, as explained in Sect. Ill the 



convex functions H(x) and ip (q) are Legendre transforms 
of each other, and the planar facets of tp(q) are associ- 
ated with vertices of H(x.), given by the second Eq.(36l, 
whereas vertices of ip(q) (i.e. the nodes of the triangula- 
tion of Fig. |6j where tp L makes contact with its convex 
hull) are associated with planar facets of i?(x), with a 
slope given by the first Eq. ( 36 1 . 



We show in Fig.[7]the resulting "Voronoi-like" tessella- 
tion of the Eulerian space obtained for a realization of the 
initial potential ipo at a given time, for the n = case. In 
agreement with the Lagrangian triangulations obtained 
in Fig. [6j we can see that in the case n = shock nodes 




FIG. 7: (Color online) The "Voronoi-like" tessellation of the 
Eulerian space associated with planar facets of the convex 
function H(x). It is the dual of the Lagrangian space tri- 
angulation of Fig. [6] We plot the diagrams in terms of the 
dimensionless scaling coordinates (Xi, X2) for the case n — 0. 



are isolated and in finite number per unit area, while void 
sizes are of order unity. In the case n = — 2 cells have 
such a small area (which keeps decreasing as we increase 
the numerical resolution) that vertices appear to cover 
the whole plane, suggesting that shock nodes are dense, 
as in the one-dimensional case. 

Finally we show in Fig. [8] the Lagrangian coordinate 
Qi (upper panel) and the velocity component U± (lower 
panel) over the Eulerian X-plane, for the case n = 0. As 
explained above, within each of the "Voronoi-like" cells 
shown in Fig. [7] the Lagrangian coordinate Q is constant, 
so that the surface Qi(X) shows a series of flat plateaus 
delimited by the boundaries of these cells, where finite 
jumps take place. Along the X 2 direction jumps in Q\ 
can be both positive and negative, with an even distri- 
bution as the system is statistically homogeneous and 
isotropic, but jumps are always positive along the X\ di- 
rection. This is a consequences of elementary properties 
of the Legendre transform: as can be seen from the first 
Eq.(36) and the fact that £f(x) is convex, we have the 
two relations (and similarly in higher dimensions) 



dqi 

dxi 



>0, 



dg 2 
dx 2 



>0, 



(47) 



which generalize the same one-dimensional property. 
However, as explained above, while in one dimension this 
non-crossing property is sufficient to reconstruct the di- 
rect Lagrangian map, q 1— > x, in higher dimensions this 
requires the use of the second Legendre transform (37 1. 



Note that by the same convexity argument we also have 



dqi 

Over the "Voronoi-like 



> 0, 



dx 2 
dq 2 



> 0. 



(48) 



cells, q being constant we can 
see from Eq.((TT|) that the velocity components are affine 
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Q<(.X 1 _X 2 ) 




FIG. 8: (Color online) The Lagrangian coordinate Qi(X) (up- 
per panel) and the velocity Ui (X) (lower panel) for a realiza- 
tion of the case n = 0. The coordinate Qi is monotonically 
increasing along direction X\ while Ui shows ramps of slope 
unity separated by downward jumps. 



functions of x. More precisely, we have for instance 
Mi(x) = (xi — qi)/t, so that within each cell u\ is con- 
stant along the x 2 direction while it grows with a slope 
1/t along the x\ direction, as can be checked in the lower 
panel of Fig. [8j At the boundaries of the cells, ui shows 
positive and negative jumps, with an even distribution, 
along the x 2 direction, and only negative jumps along the 
x\ direction, as can be seen from the behavior of q\ (x) . 



B. Merging and fragmentation 

As time goes on, the tessellation and its associated 
triangulation are bound to evolve. For any dimension, 
Lagrangian vertices (i.e. contact points between the lin- 
ear Lagrangian potential and its convex hull) at a time t 2 
form a subset of the vertices obtained at an earlier time 
t\ < £2 j as can be most easily seen from the paraboloid 
construction ( 12 ). However, it is clear that a naive merg- 



ing of neighboring triangles in the Lagrangian-space tri- 
angulation shown in Fig. [6] does not always produce a 
new triangulation, as four-sided polygons would usually 
appear. Therefore, new triangulations obtained at time 
t 2 > ti cannot be merely coarser partitions of the initial 
one. That is, even though all vertices seen at t 2 were 



already vertices at t\ , triangles seen at t 2 are not always 
the union of smaller triangles obtained at t±. This implies 
that a redistribution of matter throughout the triangu- 
lation must take place. That can only be so through 
fragmentations of shock nodes. 

In order to show how it works let us simply exam- 
ine, as in Fig. |9j successive snapshots of the triangula- 
tions/tessellations. Those presented have been obtained 
for one realization of a field of index n = 0.5 (this value 
was chosen so that cells are of similar sizes to make the 
figure easier to read but the qualitative behavior is iden- 
tical to the n — case). We display in the upper row the 
Lagrangian space triangulations and in the lower row the 
Eulerian space "Voronoi-like" diagrams, obtained at four 
successive time (from left to right). The time steps have 
been chosen so that only a few rearrangements can be 
seen from one snapshot to the other. For each output, 
we color the triangles (in Lagrangian space), and the as- 
sociated shock nodes (shown by triangular symbols in the 
lower panel in Eulerian space), that are going to be af- 
fected. Note that the relative orders of the Lagrangian 
triangles and of the Eulerian nodes match, in agreement 
with the monotonicity relations ( 47 ) and ( 48 ) . 



On these successive snapshots, only two types of tri- 
angle rearrangement can be observed: 

• a flip (2 — > 2): diagonal inversion in a quadrangle 
formed by two adjacent triangles; 

• a 3-merging (3 — > 1): merging of three adjacent 
triangles into a single one. 

These properties can be understood from the convex 
hull construction. For a generic time t, the convex hull 
is entirely made of triangular facets. This means that no 
four summits of the convex hull are coplanar. As time is 
changing, critical time values can be reached where four 
points become planar. In between two critical time val- 
ues the Lagrangian triangulation simply does not evolve 
while the associated Eulerian tessellation evolves continu- 
ously. At critical times, there are two possibilities: either 
one of the four points is within the triangle formed by the 
three others or not. The latter case corresponds to a flip 
transition; the former to a 3-merging. Here we must note 
that the existence of "flip events" , that is 2 — > 2 collisions 
with a redistribution of matter, was already noticed in [5] 
(the two-clump collision accompanied by mass exchange 
shown in their Fig. 6.21). 

Other transitions corresponding to a higher-order de- 
generacy (5 or more coplanar points) are possible in prin- 
ciple but have a vanishing probability to occur. 

Thus, we now clearly see how the system evolves with 
time through a succession of two-body and three-body 
collisions. Two-body collisions of shock nodes only re- 
distribute matter over two new shock nodes and redefine 
the Lagrangian triangulation, which allows the forma- 
tion of specific central configurations where the union of 
three neighboring triangles forms a larger triangle. Next, 
such configurations allow the removal of the central La- 
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FIG. 9: (Color online) The evolution with time of the Lagrangian (upper row) and Eulerian (lower row) tessellations, for one 
field realization with n = 0.5. We show snapshots obtained at four successive times to < t\ < t% < tz, from left to right. 

• At time to (first column) we can see in the upper panel three colored triangles which happen to form a larger triangle. They 
correspond to the three shock nodes around {X\ ~ 1.2, X2 — 1) in the lower panel, which form a triangular "Voronoi" cell. 
As seen in the second column, by time t\ these three Lagrangian triangles have merged to leave the unique larger triangle. In 
doing so the central vertex q c has been removed, while in Eulerian space the three shock nodes have merged. 

• At time t\ we have two colored triangles associated with the two shock nodes around (X\ ~ 2,X2 — 1.6) that are moving 
closer (compare with their position at to). Between time t\ and ti these two shock nodes collide and form two new shocks 
moving outward in a roughly orthogonal direction, as seen in the third column, while in Lagrangian space there has been a flip. 

• A similar process occurs between times ti and tz, with the collision and "bouncing" of the shock nodes around (Xi ~ 
2.2,A 2 ~ 1.8). 

• The new triangulation obtained between times t\ and t2 has produced the configuration with the three colored triangles 
shown at time tz, which now form a unique larger triangle (that was not the case at the earliest time to). They are associated 
with the three shock nodes around (X\ ~ 2,X2 — 1.5) which are moving closer. Thus, we recover a central configuration 
similar to the one shown at time to, and these three Lagrangian triangles and Eulerian shock nodes merge at a later time £4 
(not shown). 



grangian vertex, through the merging of the three tri- 
angles, which corresponds to a simultaneous merging of 
three shock nodes in Eulerian space. 

To conclude this paragraph we present in Fig. [10] a re- 
alization of a small fraction of a merging-fragmentation 
graph obtained for the case n = 0. The horizontal plane 
corresponds to the Eulerian x-plane while the vertical 
axis is the time (units are arbitrary). Tubes correspond 
to trajectories of halo nodes and their sections are pro- 
portional to their mass. One can see that, similarly to 
the ID case, the trajectories are straight lines in between 
the collisions taking place at critical times. In the d = 2 
case, collisions lead either to2f> 2 scatterings with mass 
exchange or to 3 — > 1 mergings. The former process is 
unknown in d = 1 case. The net result of these effects [45] 
is to progressively diminish the number density of halos 
- and to augment their average mass - in agreement with 
the scaling law 



C. Velocities in fragmentation events 

Although it was already noticed in [8], the existence 
of flip-like transitions, that is 2 — » 2 collisions with mass 
exchange in 2D, is not widely known in the cosmologi- 
cal community. In particular, despite the "geometrical 
model" , which gives rise to such events, has been studied 
in numerical works [TOl [H] in the context of the forma- 
tion of large-scale structures in cosmology, other workers 
in this broader field (who have not necessarily studied the 
problematics associated with the Burgers equation) are 
not always aware of the fact that this peculiar model im- 
plies such splitting events in dimensions two and higher 
(at least, the authors of this paper were not aware of this 
behavior before working on the present study). 

Let us point out that the Eulerian-space and 
Lagrangian-space tessellations have a different status in 
this regard. Indeed, the Voronoi-like tessellation itself, 
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FIG. 10: (Color online) A merging-fragmentation tree ob- 
tained for one realization for the case n = 0, d — 2. The bot- 
tom plane corresponds to the initial time where the Voronoi- 
like structure showing the location of the halo nodes is shown. 
The top plane corresponds to the final time. The tubes in be- 
tween show the halo trajectories. The section of each tube is 
proportional to its corresponding halo mass. One can observe 
on this picture examples of 2 o 2 halo diffusions and 3—^1 
mergings. Units are arbitrary. 



geometrical constraint (requiring that the mass distribu- 
tion is always denned by the dual Eulerian Voronoi-like 
tcsscllation/Lagrangian triangulation) implies in dimen- 
sions 2 and higher that collisions can redistribute mat- 
ter among mass clusters. In contrast, in the "standard" 
model, where mass clusters do not split, the tessellations 
are not sufficient to recover the matter distribution since 
clusters are not necessarily located at the summits of the 
Voronoi-like cells (then, one needs to integrate the conti- 
nuity equation over previous times to solve the problem) . 

In order to illustrate more clearly how the 2 — > 2 col- 
lisions proceed from a dynamical point of view, we show 
a simple example of such a flip-like transition in Fig. |11| 
Such an event was also described in Fig. 6.21 of Rcf. [8 , 
but we complete the picture by plotting the associated 
Eulerian velocity field u(x, i) and discussing in somewhat 
more details how this proceeds. This example also allows 
us to compare in Sect. |VD| below the behavior of the 
"geometrical model" with that of the "standard" model, 
using the results obtained recently by [THJ [T5] . We con- 
sider the simplest case where the initial velocity potential 
VJo(q) shows four spikes, labeled as points A, B, C and D, 
of coordinates 

L > £ > : A = (-£,0), B = (£,0), 

C=(Q,L), D = (0,-L), (49) 

with 



shown in the bottom row of Fig. [9] and in Fig. [7J is an 
Eulerian-space construction that describes the regular re- 
gions. It is thus entirely defined by the Burgers equa- 
tion in its inviscid limit, through the Hopf-Cole solution 
( fTo| . This is not so for the triangulation that lives in La- 
grangian space, shown in the upper row of Fig. [9] and m 
Fig. [6j It requires a prescription on where the matter is 
actually going, that is, one must add to the Burgers equa- 
tion a second equation for the evolution of the density 
field. One can choose the "standard" continuity equa- 
tion, or the modified continuity equation (31) associated 
with the "geometrical model" described in Sect. |III B[ 
which we study in this paper. Different models lead to 
different Lagrangian-space tessellations, since mass clus- 
ters may or may not split depending on the chosen pre- 
scription. This also means that, even though they do 
not change the Voronoi-like diagrams, that is, the cellu- 
lar structure built by the shock manifold, these different 
models put matter at different positions on this shock 
manifold. As we have already mentioned, in the "stan- 
dard" model mass clusters cannot split but they may 
leave shock nodes, as found in [THl [T5], whereas in the 
"geometrical model" mass clusters can split but they are 
always located on shock nodes (and each node is associ- 
ated with a mass cluster). Of course, the latter property 
is due to the geometrical construction of this model, as- 
sociated with the Legendre transformations and convex 
hull constructions described in Sect. |III B| As explained 
in the previous section, and illustrated in Fig. [9j this 



MA) = MB) = 0, MC) = MD) = M > 0. (50) 

Other points in the q-plane have much smaller values of 
■00 so that they do not play any role at the time of interest 
(we can take V'o(q) = ~oo everywhere outside of points 
{^4, B, C, D}). This simple case obeys at all times the two 
parity symmetries x\ <H- —x\ and xi <-> — xi. It is easiest 
to analyze the system with the help of the paraboloid 
construction (12 1. At early time (left panel in Fig. [TT|), 



paraboloids have a large curvature (i.e. they are highly 
peaked) so that Eulerian positions x in the neighborhood 
of each summit {A, B, C, D} are governed by the closest 
of these four points (i.e. paraboloids of center x make 
first contact with the closest peak among {^4, B, C, D}). 
In particular, from Eq.(ll) the velocity field u(x) out- 



flows from each of these points, as seen in left panel of 



Fig. 11 (blue arrows). Next, the "domain of influence" of 
summit A in Eulerian space (i.e. its "Voronoi-like" cell) 
is delimited by straight lines defined as the set of points x 
such that the first-contact paraboloid of center x simul- 
taneously touches summits A and B, A and C, or A and 
D (and similarly for other summits). Of course, by sym- 
metry the frontier between the A-cell and the B-ce\l is 
orthogonal to the segment (AB). These cells are shown 
by the black lines in Fig. [TT] Thus, at early times we 
have the configuration displayed in the left panel, with 
two shock nodes, shown by the big black dots, at the ver- 
tices of these Eulerian cells. Then, the matter that has 
fallen into the upper shock node comes from the upper 
triangle (ACB), whereas the lower shock node contains 
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FIG. 11: (Color online) A collision of two shock nodes that gives rise to two new shock nodes with a redistribution of mass. 



the matter from the lower triangle (ADB). This gives 
the triangulation of the Lagrangian q-space, which we 
show by the red dashed lines in the figure. Note that in 
each panel we superpose the Eulerian space (the velocity 
field marked by arrows, the "Voronoi-like" cells marked 
by solid lines and the shock nodes marked by big dots) 
and the Lagrangian space (the triangulation marked by 
the dashed lines). From these partitions of the x— and 
q— spaces we can also read the convex functions -ff(x) 
and y(q) (i.e. this gives the projection of their planar 
facets). 

As time grows the two shock nodes move closer to- 
wards the center of the figure, until the merge at time 
shown in the middle panel. At this time, the paraboloid 
of center (0, 0) makes simultaneous contact with the four 
summits {^4, B, C, D}, and the two triangles obtained at 
earlier times in the left panel merge to form a losange. 

Next, at later times we obtain the configuration shown 
in the right panel, such that the paraboloid of center 
(0, 0) only makes contact with summits C and D. Wc 
now have two shock nodes moving outward from the cen- 
ter along the horizontal axis. The Lagrangian triangle 
associated with the left (resp. right) shock node is now 
(ACD) (resp. (BCD)). Therefore, the matter within the 
losange (ACBD) has been redistributed. In particular, 
the triangle (ACB) of the left panel has been split into 
two parts, the left half going into the left shock and the 
right half going into the right shock. Thus, the unique 
central shock node obtained at time has fragmented 
into two new shock nodes and some particles that had 
coalesced at earlier times t < have separated and now 
belong to two different objects. 

It is interesting to note that it had been noticed in 
numerical simulations |14j that massive clumps seen in 
Af-body simulations of the gravitational dynamics gener- 
ally are associated with several nodes in the correspond- 
ing adhesion-model simulations (i.e. defined by the same 
initial conditions). This was interpreted as an inaccu- 
rate description of the merging process, as the Burgers 
dynamics does not take into account local gravitational 
forces. Our results show that such differences may also 
be due to the specific fragmentation process associated 



with this Burgers dynamics. 



D. Comparison with the "standard" model 

The peculiar fragmentation process described in the 
previous section arises from our prescription ( 29 1 for the 



density field, that is, from the modified continuity equa- 
tion (31). If we choose to set the right hand side of 



Eq.(31 ) to zero, that is, for finite v the density field obeys 



the standard continuity equation, we would obtain a dif- 
ferent result with no fragmentation in the inviscid limit as 
studied in [IHJQjj]. To illustrate this point, let us consider 
the configuration shown in Fig. [TTJ for a small non-zero 
viscosity, v > 0. From Eqs.Q and (10) we know that in 
the limit v — > + the velocity field at any time converges 
towards the inviscid solution displayed in Fig. [TT] The 
main difference is simply that at finite viscosity the dis- 
continuities along the shock lines (black solid lines in the 
figure) are smoothed over a small finite distance. In par- 
ticular, the central shock node obtained at time (mid- 
dle panel) has a small nonzero extension. Then, since 
the parity symmetries of the system are still exactly sat- 
isfied, we can see that at infinitesimally later times the 
particles located in the left part of this finite-size cloud 
experience a local velocity field that points towards the 
left, whereas particles located in the right part see a ve- 
locity field that points towards the right. However, if we 
use the standard continuity equation this does not lead 
to a splitting into two new symmetrical halos, moving 
further apart as time grows. Indeed, the differences be- 
tween the velocities u(x, t) of the left and right parts of 
the small central cloud go to zero with v (along with the 
size of the cloud) so that the small cloud stays in the 
middle of the new horizontal shock line in the inviscid 
limit. On the other hand, if the spikes of the initial po- 
tential tpQ at points {A, B,C, D} are not infinitesimally 
thin, two new horizontal nodes appear as in the figure 
and are fed by the matter flowing from the regular parts 
of these spikes (so that one obtains three mass clusters, 
a motionless central one which has stopped growing and 
two small outward-going ones which have just formed). 
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As mentioned earlier, the evolution of matter within 
shock lines was studied in [18l H5] . using the standard 
continuity equation for the evolution of the density field. 
It was found that the trajectories obtained in the in- 
viscid limit exist and are unique, so that trajectories 
that pass through a point at a given time coincide at all 
later times. Therefore, there is no fragmentation of ha- 
los in this prescription. On the other hand, clusters can 
stop growing and leave the shock nodes (while remaining 
on shock lines), in agreement with the discussion above 
where we explained how within this prescription a small 
halo would remain motionless at the center of the right 



the triangle 



panel of Fig. 11 This approach can be extended to more 
general Hamilton-Jacobi equations, with convex Hamil- 
tonians. Then, one again obtains a unique coalescing flow 

nana. 

The evolution with time of the shock lines themselves 
(bifurcations, transitions) in 2D and 3D was studied in 
details in [41j . for general initial conditions. In particu- 
lar, both the 2 — > 2 "flip" and 3 — > 1 merging events 
described above in the "late-time" regime for the 2D 
case correspond to the A\ space-time event of Fig. 2 in 
|41j . where the first-contact paraboloid makes simulta- 
neous contact with four "non-degenerate" points (in this 
classification a contact point is called "non-degenerate" 
if ("P XjC — tpo) has a non-zero second derivative there - in 
our case the second derivative is actually infinite). 



E. Momentum exchange through shock node 
collisions 



In the regime studied in this article, associated with 
the power-law initial conditions ([5| or the late-time stage, 
not only is all the matter content distributed within shock 
nodes but for d > 2 shocks form a unique connected set 
that spans the entire space. As a consequence, conser- 
vation of the "effective momentum" defined in Eq.(38) 



during node collisions cannot be inferred from the gen- 
eral result described in sect. |III C[ no volume V x with a 
regular boundary can be drawn around a finite number of 
nodes. We explore here how the actual momentum of the 
particles can be exchanged through shock node collisions. 

First, let us express the momentum of a given shock 
node. We note {A, B, C} the summits of its associ- 
ated triangle in Lagrangian space, which corresponds to 
a triangular facet of the convex hull ip(q) of the La- 
grangian potential. We consider time values between 
critical times, that is, a time interval during which the 
triangle {A,B,C} is left unchanged. 

The relation ( 44 1 can be written for A and B and for 
A and C . Differentiating those equations with respect 
to time t, we get the node velocity v along the vectors 
<lAB = <lA and q AC = q c - q^, 

v-q AB = ^o(cu) - V>o(qs); (51) 
v-q A c = Mqa) ~ V'o(qc)- (52) 
Furthermore the mass of the node is given by the area of 



Po , | 
m = ylQAB x q AC \ . 



(53) 



Introducing the momentum p = mv, those results can 
be recapped in a single expression. Defining the points 
{A, B, C} in a fictitious 3D space 



we obtain 



A = {qiA, q2A, Vo(qu)},---, 



V = {p, m} = y AB x AC, 



(54) 
(55) 



where AB = B — A, provided {q_AB,q_Ac} is positively 
oriented. 

Note that we use the letter "v" for the shock node 
velocity to distinguish this quantity from the Eulerian 
velocity u(x), which is discontinuous along shock lines. 
Moreover, as seen in the previous section and in Fig. 11 
in dimensions d > 2 the former is not given by the mean 
of "left" and "right" Eulerian velocities (at t > t„ the two 
shock nodes shown in the right panel of Fig. [TT] have con- 
stant finite outward velocities v whereas the horizontal 
component of u is linear over x\ and goes to zero at the 
center of the figure). Thus, we consider in this section 
the standard momentum of the particles, rather than the 
"effective momentum" introduced in Sect. MI Cl 

Let us now consider a 2 — > 2 shock collision, such as 
the one shown in Fig. [TTJ Without loss of generality, we 
can still choose the common edge (AB) before collision 
on the horizontal axis, C in the upper half plane and 
D in the lower half plane in Lagrangian space, but with 
otherwise arbitrary coordinates (i.e. the triangles need 
not be symmetrical). Then, before collision the total 3- 
momentum reads as 

V = {pabc + Pads, m A BC + m ADB } 

= V A bc + Vadb, (56) 

where we introduced the 3-momenta of the two shock 
nodes, associated with the Lagrangian triangles (ABC) 
and (ADB). At the time of the collision, paying attention 
to the orientation of the Lagrangian-space triangles, we 
obtain from Eq.(55l, 



V = 



Pa 
2 

po 
2 

Pa 
2 



ABxAC + ADx AB 
AB x VC 
VC x VA + VB x VC 



— Vdca + Vdbc- 



(57) 



At collision time, we see clearly that both the mass and 
the 2-momcntum p are conserved by the collision, as 
the 3-momenta associated with the two new shock nodes 
again sum up to V. Note that the fact that both tri- 
angles (DC A) and (DBC) are counterclockwise (so that 
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Eq.(55 1 applies with no further negative sign) comes from 
the constraint that we have a 2 — ¥ 2 collision rather than 
a 3 — > 1 event (i.e. the 2D segment [CD] intersects the 
segment L4£?]). 

Finally, let us consider a 3 — > 1 merging event, such as 
the one shown in left panel of Fig. [9] Thus, we choose a 
direct triangle (BCD) with an inner summit A. Before 
merging the 3-momenta associated with the three shock 
nodes (and the three Lagrangian-space triangles) read as 



Vabc = ^AB x AC, Vacd = ^AC x AV, 



V A 



DB 



Pa 



AV x AB, 



(58) 



whence after straightforward manipulations 



Vabc + Vacd + Vadb 



P -^BC x BV 



V 



BCD- 



(59) 



Therefore, momentum and mass are again conserved by 
the 3 — > 1 collision. 

The analysis described above, based on expression 
( 55 1 , shows that momentum is conserved in the "late- 
time regime" by shock node collisions. Note that this 
only holds at a global level: the total momentum of the 
two nodes is conserved in a fragmentation event but the 
momentum of a given node is not necessarily equal to 
the total initial momentum of the matter particles it con- 
tains. Thus, in-between arbitrary times t\ < t2, momen- 
tum is only conserved over disjoint sets of particles, such 
that particles in a given set have only been in contact 
during this time interval with particles of the same set 
(so that there has been no exchange of momentum be- 
tween different groups). 

In the d = 1 case though, there are only merging events 
and momentum conservation is therefore ensured at the 
local level, in agreement with well-known general results 

It is interesting to note that the momentum (|55| co- 



incides with Eq.(41) if the initial velocity potential ipo is 



affine over the three edges (AB), (BC) and (CA) of the 
Lagrangian-space triangle {A, B, C}. In particular, this 
means that the momentum of a given node is equal to the 
total initial momentum of the matter particles it contains 
if the initial velocity potential tpo(q) is affine over this 
triangle {A, B, C} (this also explicitly shows that both 
quantities are generically different for arbitrary potential 
tjjQ ). Then, for the power-law initial conditions ^ where 
the transition to the "late-time" regime takes place at 
an infinitesimally small time, — » + , with a charac- 
teristic scale £(£*) of the Lagrangian-space triangulation 
that goes to zero, we can see that the initial momen- 
tum is conserved (in a global sense) since the piecewise 
affine approximation of ipg defined by this triangulation 
converges to 4>o- 



VI. THREE-DIMENSIONAL DYNAMICS AND 
BEYOND 

We now briefly discuss how the results illustrated in 
the previous sections in dimensions d = 1 and d = 2 
generalize to higher dimensions, and in particular to d = 
3. 

As noticed in [8], in three dimensions the Lagrangian- 
space partition is made of tetrahedra (instead of trian- 
gles in d = 2), while the Eulerian-space partition is made 
of Voronoi cells with an arbitrary number of summits. 
Then, the analysis developed in the previous sections 
shows that the "flip" and "three-merging" events dis- 



cussed in the two-dimensional case in Sect. VB generalize 
as the following rearrangements, 

• 2 — > 3: two adjacent tetrahedra are reorganized 
into three new tetrahedra; 

• 3 — s- 2: three adjacent tetrahedra are reorganized 
into two new tetrahedra; 

• 4 — > 1: merging of four adjacent tetrahedra into a 
single one. 

These events naturally involve five vertices since, as 
seen from the paraboloid construction (12), an Eulerian- 



space shock node is associated with four Lagrangian- 
space vertices (where the paraboloid makes simultane- 
ous first-contact with ipofa)), and rearrangements occur 
when a fifth summit encounters the hyperplane defined 
by these four points. Higher-order events have a vanish- 
ingly small probability. 

The last event, 4 — > 1, which corresponds to the 
3— merging event observed in d — 2, occurs when a 
Lagrangian-space summit which is located within the 
tetrahedron built by four neighboring summits is re- 
moved from the convex hull, as in the first step shown in 
left panel in Fig. [9] Indeed, before its removal, this inte- 
rior summit leads to a splitting of the larger tetrahedron 
into four distinct tetrahedra, associated with four shock 
nodes. After the removal, only the embedding tetrahe- 
dron is left, which corresponds to a single shock node that 
contains all the mass associated with the four previous 
shock nodes. 

When the convex hull of the five Lagrangian vertices 
is not a tetrahedron (i.e., no point is located within the 
tetrahedron formed by the other four summits) , one can- 
not build a single tetrahedron by complete merging, so 
that, as in the two-dimensional events shown in middle 
panels of Fig. |9j one can only observe rearrangements of 
the matter distribution into new tetrahedra. This leads 
to both events 2 — > 3 and 3 — > 2. 

We illustrate the first case, 2 — > 3, in Fig. [l2j This 
corresponds to Fig. [TT] associated with 2 — > 2 events in 
d = 2. We use the same notations and a similar symmet- 
rical configuration. Thus, we have three Lagrangian sum- 
mits in the equatorial plane, z = 0, which form an equi- 
lateral triangle of side I with the same value ipo qU y an d 
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t <t, t >t. 




FIG. 12: (Color online) A collision of two shock nodes that gives rise to three new shock nodes, with a redistribution of mass, 
in the three dimensional case. Left panel: at early times we have two shock nodes (black dots), moving closer to the equatorial 
plane, which contain the mass associated with the upper and lower tetrahedra (red dashed lines). Right panel: after collision 
at t* we have three shock nodes moving apart in the equatorial plane, while in Lagrangian space the triangular bipyramid built 
by the two previous tetrahedra has been split into three new tetrahedra, with a common edge set by the segment joining the 
upper an lower summits (vertical dashed line). 



two symmetric summits on the vertical axis at a larger 
distance L with the same value tpQ^ d > 4>o qn - We again 
superpose the Eulerian-space Voronoi cells (with edges 
shown by the black solid lines) and the Lagrangian-space 
tessellation built from elementary tetrahedra (dashed red 
lines). Then, at early times, using again the Hopf-Cole 
paraboloid solution ( 12 ), each Lagrangian summit is con- 



tained within its associated Voronoi cell and for L large 
enough (or at small enough t) the central point (0, 0, 0) 
"sees" the three equatorial summits. This leads to the 
partition shown in the left panel with two shock nodes. 

As time increases, since tpo^ > "0o qu the Voronoi cells 
associated with the two upper and lower points expand 
and eventually make contact at the center, at the time 
t* when the two shock nodes collide. Afterwards, these 
two Voronoi cells have a common triangular facet in the 
equatorial plane, the cells associated with the three equa- 
torial summits being pushed outward, and we have three 
outward-moving shock nodes. In Lagrangian space, the 
two symmetric upper and lower tetrahedra shown in the 
left panel, which merge into a unique volume at t = t*, 
have been split into three new tetrahedra, with a com- 
mon edge given by the segment that joins the upper and 
lower summits. 



As for the 2D case shown in Fig. 11 it is interesting 
to note that the "standard" model would give a different 
behavior. Indeed, extending the discussion of Sect. |VD| 
we can see that using the standard continuity equation 
we would obtain a single mass cluster that remains mo- 
tionless in the center of the figure. Thus, there is no frag- 
mentation but the cluster leaves the shock nodes while 
staying within the shock manifold. Note that this mass 
cluster is no longer located on shock lines either, since it 
sits at the center of the horizontal triangular facet of the 
shock manifold. This explicitly shows that within this 
"standard" model matter can be distributed anywhere 



on the shock manifold, that is, not only on shock nodes 
or shock lines, and that one needs to know the evolution 
of the system over all previous times. 

It is clear that the reverse event, 3 — >• 2, can be ob- 
tained in a similar fashion, by choosing the upper and 
lower summits close to the equatorial plane (L < i) with 

€ /d < C U - 

These processes are naturally expected to generalize 
to higher dimensions, where the Lagrangian-space trian- 
gulation is built from (d + l)-summit cells. Then, com- 
plete mergings are associated with (d+1) — > 1 collisions 
in Eulerian space, while lower-order collisions, 2 — > d, 
3 (d — 1), d 2, are associated with Lagrangian- 
space rearrangements. The latter lead to a redistribution 
of matter, so that particles which had coalesced into the 
same shock node at an earlier time can be separated into 
distinct objects. Moreover, such collisions also change 
the number of mass clusters in a very specific manner, 
according to these rules. 



VII. CONCLUSIONS, DISCUSSION 

The Burgers equation, and its Hopf-Cole solution, only 
determines the evolution of the velocity field and to cou- 
ple this to a transportation of matter one must add an 
equation for the evolution of the density field. A nat- 
ural choice would be to use the "standard" continuity 
equation, as in [T5J [TH] , but it should then be integrated 
numerically over time. In the inviscid limit, it also leads 
to behaviors that, at least in a cosmological context, are 
not necessarily realistic. 

Here we rather explore an alternative choice, which we 
call the "geometrical model" , that fully takes advantage 
of the Hopf-Cole solution to define the matter distribu- 
tion from a geometrical construction. It is based on Leg- 
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endre transformations and convex hull constructions, so 
that the matter distribution is associated with dual Eu- 
lerian and Lagrangian space tessellations [HI HO] • In reg- 
ular regions and in the inviscid limit, both approaches 
coincide. The geometrical model however corresponds to 
a non-standard continuity equation for the density field 
that affects the mass behavior within the shock mani- 
folds. 

The peculiarities of this model are best revealed in the 
late time behavior of the density field or for power-law 
spectra initial conditions. Then, in the inviscid limit that 
we investigate in this article, the matter is entirely con- 
tained in shock nodes, which have gathered the matter 
that was initially in segments, triangles, or tetrahedra, 
for respectively the d = 1, d = 2, or d = 3 cases. In 
any dimension, those objects form a partition of the La- 
grangian space. 

The most striking result of these investigations is the 
mechanism with which the matter is rearranged in ha- 
los of growing mass as the system is evolving with time. 
Indeed, for dimensions greater or equal to 2, these rear- 
rangements follow a complex pattern of successive "flips" 
and "mergings". In d = 2, we explicitly show that these 
"flips" correspond in Lagrangian space to the rearrange- 
ments of two triangles into two new triangles, and in Eu- 
lerian space to a collision of two shock nodes giving rise 
to two new outward- moving shock nodes. The "merging" 
events correspond in Lagrangian space to the merging of 
three triangles into a single larger triangle, and in Eu- 
lerian space to the simultaneous merging of three halos 
into a single node. More complex patterns are expected 
for higher dimensional cases. 

We have also described how the "flips" or "fragmenta- 
tion" events, associated with collisions that give rise to 
several new shock nodes, lead to a redistribution of mo- 
mentum over outgoing nodes. This implies that momen- 
tum is only conserved in a global sense in d > 2 (whereas 
conservation of momentum also holds in a stronger local 
sense in d = 1, where there are only merging events). 
Thus, in this regime the inviscid limit of the Burgers 



equation leads to a specific set of collision rules (see note 
[4"B]). such as {2 ->• 2,3 -> 1} in 2D and {2 -> 3, 3 -> 
2,4 — > 1} in 3D. In the present case, even though shock 
nodes are infinitesimally thin, we have seen that one must 
take into account up to (d + l)-body collisions in dimen- 
sion d, which still occur at a finite rate. 

This peculiar behavior is due to the geometrical con- 
struction that defines the matter distribution and that 
allows to integrate the equations of motion for both the 
velocity field (using the standard Hopf-Cole solution of 
the Burgers equation) and the density field (by definition 
of this geometrical construction) . This is evidently a very 
convenient property that has motivated this study. How- 
ever, as shown here, this leads to a singular inviscid limit, 
in the sense that particles that share the same location 
at a given time can separate at later times - such singu- 
lar behaviors are associated with collisions between shock 
nodes - so that it does not lead to a genuine "adhesion 
model". This result emphasizes the fact that the trans- 
portation of matter which can be associated to the Burg- 
ers equation is a non-trivial process, whether one uses 
the "standard" or the "geometrical" model, and that in 
the latter case one must be aware of the peculiar merging 
and fragmentation events that take place as shock nodes 
collide and reorganize the shock manifold. 
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simple example where these conditions are identically 
satisfied is the case where the initial velocity potential 

is separable, V'o(q) = E« (?*)> anc ^ the dynamics can 
be factorized in terms of d one-dimensional dynamics. 
In generic systems they are only satisfied at peculiar 
isolated points (note that to obey the last condition it is 
not sufficient to have a parity sym metr y Xi —Xi). 
[44] Rather than using the prescription (361, some numerical 
simulations try to infer the distribution of matter directly 
from the Hopf-Cole solution. For instance, [131 115] first 
compute the velocity field over a grid of times from the 
explicit solution (TsJ) , with a small finite viscosity v, and 
second integrate the particle orbits along this velocity 
field . Alternatively, working directly in the in visc id limit 
( |10[ ), [121114] use the paraboloid construction ( 12 l to con- 
struct the "skeleton" of the Eulerian-space density field 
(i.e. the Voronoi cells discussed in sect. |V A] in this paper) . 
Thus, in 2D, they classify the Eulerian-space positions x 
(on a grid) in three classes, (1) regular points where there 
is only one first-contact point q between the paraboloid 
and the initial potential ipo, (2) singular points with two 
simultaneous first-contact points, and (3) singular points 
with three first-contact points. Points (2) make straight 
segments, connected at vertices (3), that build a Voronoi- 
like tessellation with cells that contain the regular points 
(1). Then, from class (1) one identifies regular regions in 
Lagrangian space, the complement of which defining the 
singular regions which have fallen into shocks (i.e. the 
filaments and nodes of the "skeleton"). Next, shocked 
matter is distributed on the skeleton by following the 
motion in the regular regions (blue arrows in Fig. 
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this paper). This actually corresponds to an indirect im- 
plementation of the prescription (36 1 used in this article. 



symmetric it can be diagonalized by an orthogo- 
nal matrix. In this orthonormal basis, denoting Xk 



[45] It can also be remarked that critical time values are 
strongly correlated, as in the center of the figure few 
subsequent fragmentations/mergings happen on a rather 
short time scale as compared to the typical time scale of 
this figure. We did not try to investigate in more details 
the time correlations of these critical events. 

[46] They are at variance with familiar discrete systems, 
where in the limit of small particle radii the dynamics is 
usually dominated by two-body collisions as higher-order 
collisions are vanishingly rare. 



